rm(list = ls())
#作业1
#导入并处理数据
m <- read.csv("./ad.csv")
rownames(m) <- m$ID
m <- m[,-1:-2]
m[m==0] <- NA
#创建数据框储存结果
p_value <- data.frame(p=1:100)
rownames(p_value)<- rownames(m)
#对数化
x1 <- log2(m)
#提取分组信息colnames(x1) <- group
group <- gsub("\\d","",colnames(x1))
table(group)
#循环
for (i in 1:nrow(x1) ){
  proi<- x1[i,]
  proi.1 <- unlist(proi)
#判断各组数据是否有至少3个样本是有数值的
ctl <- proi.1[group=="ctl"]
ad <- proi.1[group=="ad"]
asym <- proi.1[group=="asym"]
ctlnum = sum( !is.na(ctl))
adnum = sum(!is.na(ad))
asymnum = sum(!is.na(asym))
#对于判定结果进行方差分析或者直接输出NA
if(asymnum >3 & adnum>3 & ctlnum >3){
  anovaresult <- oneway.test(proi.1~group)
  p <- anovaresult$p.value
} else p = NA
#将p值存入对应的结果框
p_value[i,1] <- p
}



#作业2
#加载包
library(clusterProfiler)
library(org.Hs.eg.db)
library(ggplot2)
#加载volcano.RData数据
load("./data/volcano.RData")
#提取p<0.05的蛋白质对应的基因列表。
genes <- prostat$ID[prostat$P < 0.05 ]
#进行GO富集分析，使用clusterProfiler的enrichGO函数。
GO <- enrichGO(
      gene = genes,        # 输入基因列表
      OrgDb = org.Hs.eg.db,     # 使用人类基因注释数据库
      keyType = "SYMBOL",     # 输入基因ID类型
      ont = "BP",               # 分析生物过程（Biological Process）
      pAdjustMethod = "BH",     # 使用Benjamini-Hochberg方法校正p值
      pvalueCutoff = 0.05,      # 显著性阈值：p值 < 0.05
      qvalueCutoff = 0.05,      # 显著性阈值：q值 < 0.05     
    )
head(GO)
#绘制气泡图
dotplot(
  GO,
  showCategory = 20,
  label_format=100,
  font.size    = 10,
  color        = "pvalue",    # 颜色映射为p值
  size         = "Count",     # 点大小映射为基因数量
 
  title        = "GO Enrichment Analysis "
) +theme(axis.text.y = element_text(margin = margin(r=5))
)
# 保存富集结果
write.csv(as.data.frame(GO), "GO_enrichment_results.csv")

# 保存气泡图
ggsave("GO气泡图.png", width = 12, height = 8, dpi = 300)
